## Function for drawing groups using blockTools and Random Allocation
## Ryan T. Moore
## 11 July 2014

draw_groups <- function(n = 24, num.groups = 6, bl.vars = "isMale", n.rands = 100){
  store.avg.bt <- store.avg.cr <- store.prop.allmenwomen <- NULL
  store.pvalFbl <- store.pvalKWbl <- store.pvalFcr <- store.pvalKWcr <- matrix(NA, n.rands, length(bl.vars))
  for(i in 1:n.rands){
    x <- data.frame(id = 1:n, isMale = sample(rep(0:1, n/2)), year = sample(1:2, n, rep = TRUE), midterm = round(rnorm(n, m = 80, sd=10)), polknowl = rbinom(n, 10, .7), distance = round(rlnorm(n, m = 5, sd = 1.2)))
    b <- block(x, id = "id", block.vars = bl.vars, n.tr = num.groups) 
    a <- assignment(b)
    x$cr <- sample(rep(1:num.groups, n/num.groups))
    x$bl <- x$bla <- NA
    for(row.idx in 1:nrow(x)){
      x$bl[row.idx] <- ceiling(which(b$blocks[[1]] == x$id[row.idx])/nrow(b$blocks[[1]]))
      x$bla[row.idx] <- ceiling(which(a$assg[[1]] == x$id[row.idx])/nrow(a$assg[[1]]))
    }
    store.gp.means.bt <- store.gp.means.cr <- NULL
    for(gp.idx in 1:num.groups){
      store.gp.means.bt <- append(store.gp.means.bt, mean(x[x$id %in% a$assg[[1]][, gp.idx], "isMale"]))
      store.gp.means.cr <- append(store.gp.means.cr, mean(x[x$cr == gp.idx, "isMale"]))
    }
    store.avg.bt[i] <- mean(store.gp.means.bt == .5) ## fraction evenly split, 1 dataset
    store.avg.cr[i] <- mean(store.gp.means.cr == .5) 
    store.prop.allmenwomen[i] <- mean(store.gp.means.cr %in% c(0, 1)) ## fraction of homog. groups, 1 dataset
    for(var.idx in 1:length(bl.vars)){
      form.b <- formula(paste(noquote(bl.vars[var.idx]), " ~ bla", sep = ""))
      form.c <- formula(paste(noquote(bl.vars[var.idx]), " ~ cr", sep = ""))
      store.pvalFbl[i, var.idx] <- oneway.test(form.b, data = x)$p.value
      store.pvalFcr[i, var.idx] <- oneway.test(form.c, data = x)$p.value
      store.pvalKWbl[i, var.idx] <- kruskal.test(form.b, data = x)$p.value
      store.pvalKWcr[i, var.idx] <- kruskal.test(form.c, data = x)$p.value
    }
  }
  output <- list(blocked = store.avg.bt, crand = store.avg.cr, allmw = store.prop.allmenwomen, pvFbl = store.pvalFbl, pvFcr = store.pvalFcr, pvKWbl = store.pvalKWbl, pvKWcr = store.pvalKWcr)
  output$call <- match.call()
  return(output)
}